Generic Chemistry

Generic chemical equilibrium and kinetic reactions following the standard approach used by reaction-transport codes such as PHREEQ and CrunchFlow, see eg (Steefel et al., 2015). This exploits timescale separation between "fast" (assumed instantaneous) chemical equilibrium reactions, and "slow" kinetic reactions or transport.

  • The chemical system is represented by a small number of totals

(or components) and an equal number of primary species concentrations, with secondary species concentrations calculated from the primary species via a set of equilibrium reactions. Primary species concentrations are determined by solving the set of algebraic equations given by the constraints on total concentrations.

  • Kinetic reactions (with any species, primary or secondary as reactants) are then written with totals as products.
  • Bulk transport (eg ocean advection or eddy diffusivity) transports totals. Molecular diffusivity (eg in a sediment) transports primary or secondary species and accumulates fluxes into totals.

Reservoirs

PALEOaqchem.ReservoirsAq.ReactionConstraintReservoirType
ReactionConstraintReservoir

A primary species and (algebraic) constraint on a corresponding total or component.

The primary species concentration or amount is defined as a PALEO State Variable, which depending on the primary_variable parameter, may be:

  • Primary_conc: (mol m-3)
  • Primary: (mol)
  • Primary_pconc: -log 10 (concentration (mol kg-1))

The corresponding R_constraint_conc or R_constraint (mol) defining the algebraic constraint on the corresponding total (for use by the numerical solver) is defined as a PALEO Constraint Variable.

This ReactionConstraintReservoir would usually be used in combination with a ReactionReservoir that provides the required total component concentration or amount as an ODE variable (where as usual reaction source and sink fluxes are applied to the corresponding _sms variable). Depending on the constraint_variable parameter, the total component may be supplied as either a per-cell concentration or amount:

  • R_conc: (mol m-3)
  • R: (mol)

Equilibrium reactions defining secondary species should add their contributions to the total to R_calc (mol). A primary species contribution R_calc += primary_total_stoich * Primary_conc * primary_volume is added to R_calc (where for the usual case parameter primary_total_stoich should be set to 1.0). Primary species contributions to other totals can be included by setting the primary_other_components parameter.

The numerical solver then solves for the primary species (and hence the secondary species concentrations) that (depending on the constraint_variable parameter) satisfy one of:

0 = R_constraint_conc = R_conc - R_calc/volume
0 = R_constraint = R - R_calc

Volume conversions

The total species concentration R_conc and primary species concentration Primary_conc use (potentially different) volume conversions provided in volume and primary_volume respectively.

This allows for cases eg equilibrium partitioning between solute and solid phases by surface complexation, where R_conc refers to a cell total volume, and Primary_conc to a solute concentration.

Parameters

  • primary_total_stoich[Float64]=1.0, default_value=1.0, description="stoichiometric factor R_calc_conc += primary_total_stoich * Primary_conc"
  • primary_other_components[Vector{String}]=String[], default_value=String[], description="contribution of primary species to other element or component total concentrations"
  • primary_variable[String]="concentration", default_value="concentration", allowed_values=["concentration", "amount", "p_concentration"], description="units for primary variable"
  • constraint_variable[String]="concentration", default_value="concentration", allowed_values=["concentration", "amount"], description="units for constraint variable"

Methods and Variables for default Parameters

  • do_constraintreservoir_primary
    • R_calc (mol), VT_ReactContributor, description="contributions to total R_calc_conc (NB: a total, not concentration, to generalize to multiphase eqb)"
    • primary_volume –> volume (m3), VT_ReactDependency, description="cell volume (as used by Primary_conc)"
    • Primary_conc (mol m-3), VT_ReactDependency, VF_State, description="concentration of primary species"
  • do_constraintreservoir_constraint
    • R_calc (mol), VT_ReactTarget, description="contributions to total R_calc_conc (NB: a total, not concentration, to generalize to multiphase eqb)"
    • R_constraint_conc (mol m-3), VT_ReactContributor, VF_Constraint, description="algebraic constraint on R_conc (= 0)"
    • R_conc (mol m-3), VT_ReactDependency, description="total R_conc"
    • volume (m3), VT_ReactDependency, description="cell volume (as used by total variable)"
source
PALEOaqchem.ReservoirsAq.ReactionImplicitReservoirType
ReactionImplicitReservoir

A primary species and corresponding total or component as an 'implicit' ODE variable.

This provides an implementation of the 'Direct Substitution Approach' to chemical speciation, where the total or component is a function of the primary species concentration.

The primary species concentration or amount is defined as a PALEO StateTotal Variable, which depending on the primary_variable parameter, may be:

  • Primary_conc: (mol m-3)
  • Primary: (mol)
  • Primary_pconc: -log 10 (concentration (mol kg-1))

The corresponding total component R_conc or R is defined as a PALEO Total Variable, which depending on the constraint_variable parameter, may be provided to the solver either as a per-cell concentration or amount:

  • R_conc = R_calc/volume: (mol m-3)
  • R = R_calc: (mol)

Equilibrium reactions defining secondary species should add their contributions to the total to R_calc (mol). A primary species contribution R_calc += primary_total_stoich * Primary_conc * primary_volume is added to R_calc (where for the usual case parameter primary_total_stoich should be set to 1.0). Primary species contributions to other totals can be included by setting the primary_other_components parameter.

Source - sink fluxes eg kinetic reactions should be added to R_sms (mol yr-1) defined as a PALEO Deriv Variable.

Volume conversions

The total species concentration R_conc and primary species concentration Primary_conc use (potentially different) volume conversions provided in volume and primary_volume respectively.

This allows for cases eg equilibrium partitioning between solute and solid phases by surface complexation, where R_conc refers to a cell total volume, and Primary_conc to a solute concentration.

Parameters

  • primary_total_stoich[Float64]=1.0, default_value=1.0, description="stoichiometric factor R_calc_conc += primary_total_stoich * Primary_conc"
  • primary_other_components[Vector{String}]=String[], default_value=String[], description="contribution of primary species to other element or component total concentrations"
  • primary_variable[String]="concentration", default_value="concentration", allowed_values=["concentration", "amount", "pconcentration"], description="units for primary variable (specifies Primary\conc, Primary, Primary_pconc as StateTotal variable)"
  • total_variable[String]="concentration", default_value="concentration", allowed_values=["concentration", "amount"], description="units for total variable (specifies R_conc, R as Total variable)"
  • total[Bool]=false, default_value=false, description="true to calculate R_total = sum(R)"

Methods and Variables for default Parameters

  • do_implicitreservoir_primary
    • R_calc (mol), VT_ReactContributor, description="contributions to total R_calc_conc (NB: a total, not concentration, to generalize to multiphase eqb)"
    • primary_volume –> volume (m3), VT_ReactDependency, description="cell volume (as used by Primary_conc)"
    • Primary_conc (mol m-3), VT_ReactDependency, VF_StateTotal, description="concentration of primary species"
  • do_implicitreservoir_sms
    • R_sms (mol yr-1), VT_ReactTarget, description="total or component R source - sink"
    • R_conc_sms (mol m-3 yr-1), VT_ReactContributor, VF_Deriv, description="total or component R_conc source - sink"
    • volume (m3), VT_ReactDependency, description="cell volume (as used by total variable)"
  • do_implicitreservoir_total
    • R_calc (mol), VT_ReactTarget, description="contributions to total R_calc_conc (NB: a total, not concentration, to generalize to multiphase eqb)"
    • volume (m3), VT_ReactDependency, description="cell volume (as used by total variable)"
    • R (mol), VT_ReactProperty, description="total or component R"
    • R_conc (mol m-3), VT_ReactContributor, VF_Total, description="total or component R_conc"
source
PALEOaqchem.ReservoirsAq.ReactionAqConcSumType
ReactionAqConcSum

A sum of concentration variables (eg to get an element total)

Parameters

  • vars_to_add[Vector{String}]=["2\myvar", "myothervar", "-1\mythirdvar"], default_value=["2\myvar", "myothervar", "-1\mythirdvar"], description="vector of variable names to add, eg [2*myvar, myothervar, -1*mythirdvar]"
  • add_to_sum_volume[Bool]=false, default_value=false, description="true to also add to a 'sum' Variable += 'sum_conc * volume"
  • define_sum_volume[Bool]=false, default_value=false, description="only if 'add_to_sum_volume == true': true to also define the 'sum' Variable"

Methods and Variables for default Parameters

  • do_aqconcsum
    • sum_conc (mol m-3), VT_ReactProperty, description="sum of specified variables"
    • myvar (), VT_ReactDependency, description=""
    • myothervar (), VT_ReactDependency, description=""
    • mythirdvar (), VT_ReactDependency, description=""
source

Equilibrium reactions

PALEOaqchem.GenericReactions.ReactionAqEqbType
ReactionAqEqb

Define a new equilibrium species N

N + a A <--> b B + c C

[N] = K_eqb'^K_power ([B]^b [C]^c) / ([A]^a)

where to convert density units for K_eqb:

K_eqb' = K_eqb * rho_ref^K_density_power

The first name in the Reactants list is the new species concentration N: other species concentrations in Reactants and Products lists must already be defined elsewhere in the model configuration.

The contribution of the new species to element totals or components is defined by the N_components vector, which may be empty eg to just calculate an Omega or a gas partial pressure etc.

Examples

Gas partial pressure from concentration

solubility_H2:
    class: ReactionAqEqb
    parameters:
        Reactants:  ["pH2"]
        Products:  ["H2_conc"]
        K_eqb:      7.8e-1   # mol m-3 atm-1 at 298.15 K (Henry's law coefficent)
        K_power:   -1.0  #  pH2 = H2_conc * K_eqb^-1
    variable_attributes:
        pH2%units:   atm

Compilation of Henry's law coefficients: https://www.henrys-law.org/ which is R. Sander: Compilation of Henry's law constants (version 5.0.0) for water as solvent, Atmos. Chem. Phys., 23, 10901-12440 (2023), doi:10.5194/acp-23-10901-2023

Unit conversions: 1 mol m^-3 Pa-1 = 1.01325e5 mol m-3 atm-1

Parameters

  • Reactants[Vector{String}]=["N\conc", "A\conc^2"], default_value=["N\conc", "A\conc^2"], description="concentrations or activities of new species followed by other reactants, write powers as X^2 etc"
  • Products[Vector{String}]=["B\conc", "C\conc"], default_value=["B\conc", "C\conc"], description="concentrations or activities of products, write powers as X^2 etc"
  • K_eqb[Float64]=0.00018629779999999998, default_value=0.00018629779999999998, description="equilibrium constant"
  • K_density_power[Float64]=0.0, default_value=0.0, description="multiple K_eqb * rho_ref^K_density_power to convert units: 0.0 for K_eqb in mol m-3, 1.0 for K_eqb in mol kg-1, etc"
  • K_power[Float64]=-1.0, default_value=-1.0, description="exponent of K_eqb"
  • N_components[Vector{String}]=String[], default_value=String[], description="contribution of new species to element or component total eg '["2*TN_calc"]' to add 2\*N\_conc\*volume to TN\_calc (or empty vector to just define an Omega)"

Methods and Variables for default Parameters

  • do_aqeqb
    • N_conc (mol m-3), VT_ReactProperty, description="aqueous concentration or activity"
    • A_conc (mol m-3), VT_ReactDependency, description="aqueous concentration or activity"
    • B_conc (mol m-3), VT_ReactDependency, description="aqueous concentration or activity"
    • C_conc (mol m-3), VT_ReactDependency, description="aqueous concentration or activity"
source

Kinetic reactions

PALEOaqchem.GenericReactions.ReactionAqKineticType
ReactionAqKinetic

Define a kinetic reaction with rate dependent on concentrations

a A + b B --> c C + d D

Rate (default) is:

R = K * [A] * [B]

where this can be modified to different functional form by defining a vector of Rate_functions to apply to each concentration variable.

Parameters Reactants and Products should be the vectors of stoichiometry * <name> of (total) species to accumulate fluxes into <name>_sms variables.

Parameter Reactant_concs should be an empty vector to take default concentration variable names from Reactants, or a Vector of names to specify concentration species names explicitly (required when eg where Reactants refers to totals which are partioned into multiple species).

Parameters

  • Reactants[Vector{String}]=["A", "2\B"], default_value=["A", "2\B"], description="reactant species"
  • Products[Vector{String}]=["2\C", "D"], default_value=["2\C", "D"], description="product species"
  • Reactant_concs[Vector{String}]=String[], default_value=String[], description="names of concentration variables to calculate rate eg '[`"A_conc"]' etc, empty vector to used defaults from 'Reactants' eg 'A_conc', 'B_conc' ..."
  • Rate_functions[Vector{String}]=String[], default_value=String[], allowed_values=["linear", "sqrt", "monod"], description="functional form for rate function of each concentration (empty vector for default 'linear')"
  • K[Float64]=NaN, default_value=NaN, description="rate constant"
  • K_lim[Float64]=NaN (mol m-3), default_value=NaN, description="limiting concentration for 'monod' rate function"

Methods and Variables for default Parameters

  • do_aqkinetic
    • redox_A_B_C_D (mol yr-1), VT_ReactProperty, description="rate variable"
    • A_conc (mol m-3), VT_ReactDependency, description="aqueous concentration or activity"
    • B_conc (mol m-3), VT_ReactDependency, description="aqueous concentration or activity"
    • volume (m3), VT_ReactDependency, description="cell solute volume"
  • RateStoich_redox_A_B_C_D
    • redox_A_B_C_D (mol yr-1), VT_ReactDependency, description="rate variable"
    • [A_sms] (mol yr-1), VT_ReactContributor, description="generated by RateStoich rate=redox_A_B_C_D"
    • [B_sms] (mol yr-1), VT_ReactContributor, description="generated by RateStoich rate=redox_A_B_C_D"
    • [C_sms] (mol yr-1), VT_ReactContributor, description="generated by RateStoich rate=redox_A_B_C_D"
    • [D_sms] (mol yr-1), VT_ReactContributor, description="generated by RateStoich rate=redox_A_B_C_D"
  • totals
    • redox_A_B_C_D (mol yr-1), VT_ReactDependency, description="rate variable"
    • redox_A_B_C_D_total (mol yr-1), VT_ReactProperty, description="total rate variable"
source

Precipitation-dissolution reactions

PALEOaqchem.GenericReactions.ReactionAqPrecipDissolType
ReactionAqPrecipDissol

Define a precipitation and dissolution reaction for solid S

a A + b B <--> s S + d D

Rate for the precipitation and dissolution reactions are:

R_precip = K_precip * (Ω - 1)               (Ω > 1)
R_dissol = K_dissol * S_conc * (1 - Ω)      (Ω < 1)

Parameters Reactants and Products should be the vectors of stoichiometry * name of (total) species to accumulate fluxes into _sms variables.

Solid_conc should be the name of the concentration variable for S, or an empty string to use default 'S_conc'.

Parameters

  • Reactants[Vector{String}]=["A", "2\B"], default_value=["A", "2\B"], description="reactant species"
  • Products[Vector{String}]=["S", "0.5\D"], default_value=["S", "0.5\D"], description="product species (solid S first)"
  • Solid_conc[String]="", default_value="", description="name of solid S concentration variable (empty string to default to 'S_conc')"
  • K_precip[Float64]=0.0 (mol m-3 yr-1), default_value=0.0, description="rate constant for precipitation reaction"
  • K_dissol[Float64]=0.0 (yr-1), default_value=0.0, description="rate constant for dissolution reaction"
  • dissol_rolloff_conc[Float64]=0.0 (mol m-3), default_value=0.0, description="limiting concentration below which dissolution rate is rolled off to zero as Solid_conc^2"

Methods and Variables for default Parameters

  • do_aqprecipdissol
    • precipdissol_A_B_S_D (mol yr-1), VT_ReactProperty, description="rate variable"
    • S_conc (mol m-3), VT_ReactDependency, description="solid concentration or activity"
    • Omega (), VT_ReactDependency, description="saturation state"
    • volume (m3), VT_ReactDependency, description="cell solid phase volume"
  • RateStoich_precipdissol_A_B_S_D
    • precipdissol_A_B_S_D (mol yr-1), VT_ReactDependency, description="rate variable"
    • [A_sms] (mol yr-1), VT_ReactContributor, description="generated by RateStoich rate=precipdissol_A_B_S_D"
    • [B_sms] (mol yr-1), VT_ReactContributor, description="generated by RateStoich rate=precipdissol_A_B_S_D"
    • [S_sms] (mol yr-1), VT_ReactContributor, description="generated by RateStoich rate=precipdissol_A_B_S_D"
    • [D_sms] (mol yr-1), VT_ReactContributor, description="generated by RateStoich rate=precipdissol_A_B_S_D"
  • totals
    • precipdissol_A_B_S_D (mol yr-1), VT_ReactDependency, description="rate variable"
    • precipdissol_A_B_S_D_total (mol yr-1), VT_ReactProperty, description="total rate variable"
source

Particulate fluxes

PALEOaqchem.Particle.ReactionParticleDecayType
ReactionParticleDecay

Decay (eg organic matter to remineralization) at rate 1.0/decay_timescale of eg an organic matter dissolved/particulate phase in Reservoir Particle, to decayflux. Particle may have isotope_type. The Reservoir Particle may have the :vsink attribute set to represent a marine sinking particulate phase.

Parameters

  • decay_timescale[Float64]=0.5 (yr), default_value=0.5, description="particle decay timescale"
  • decay_threshold[Float64]=-Inf (mol m-3), default_value=-Inf, description="particle decay concentration below which decay stops"
  • show_decayrate[Bool]=false, default_value=false, description="true to create a 'decayrate' variable to record decay rate"
  • field_data[DataType]=PALEOboxes.ScalarData, default_value=PALEOboxes.ScalarData, allowed_values=Type[PALEOboxes.ScalarData, PALEOboxes.IsotopeLinear], description="disable / enable isotopes and specify isotope type"

Methods and Variables

  • do_particle_decay
    • Particle (mol), VT_ReactDependency, description="Particle amount"
    • Particle_sms (mol yr-1), VT_ReactContributor, description="Particle source-sink"
    • decayflux (mol yr-1), VT_ReactContributor, description="Particle decay flux"
source
PALEOaqchem.Particle.ReactionFluxToComponentsType
ReactionFluxToComponents

Distribute a single input_flux (eg an organic matter flux) to output_fluxnames components with stoichiometry output_fluxstoich. input_flux may have an isotope type (set by field_data) in which case the isotopic composition will be sent to (usually one) output_fluxname with ::Isotope suffix.

Parameters

  • outputflux_prefix[String]="", default_value="", description="Prefix for output flux names"
  • outputflux_names[Vector{String}]=["Corg", "N", "P"], default_value=["Corg", "N", "P"], description="Suffixes for output flux names. Use ::Isotope suffix to identify a flux with 'isotope_type'"
  • outputflux_stoich[Vector{Float64}]=[106.0, 16.0, 1.0], default_value=[106.0, 16.0, 1.0], description="Stoichiometry for output fluxes relative to input flux"
  • allow_unlinked_outputflux[Bool]=false, default_value=false, description="true to allow (and ignore) unlinked outputflux Variables, false to error if any outputflux Variable is not linked"
  • field_data[DataType]=PALEOboxes.ScalarData, default_value=PALEOboxes.ScalarData, allowed_values=Type[PALEOboxes.ScalarData, PALEOboxes.IsotopeLinear], description="disable / enable isotopes and specify isotope type"

Methods and Variables for default Parameters

  • do_flux_to_components
    • inputflux (mol yr-1), VT_ReactTarget, description="input flux"
    • Corg (mol yr-1), VT_ReactContributor, description="flux Corg"
    • N (mol yr-1), VT_ReactContributor, description="flux N"
    • P (mol yr-1), VT_ReactContributor, description="flux P"
source

Co-precipitation

PALEOaqchem.CoPrecip.ReactionPACoPrecipType
ReactionPACoPrecip

Co-precipitation of P (eg iron-bound phosphorus) with A (eg Fe oxide) formation

P -> P=A

at a rate gamma * A_formation_rate_<n>, with limitation at low P concentration

P_components defines solution totals or components that should be modified when P is consumed (eg ["-1*P"] to remove "P" from solution)

Parameters

  • A_rate_stoich_factors[Vector{Float64}]=[1.0], default_value=[1.0], description="stoichiometry factor to multiply each A formation rate variable to convert to mol A"
  • gamma[Float64]=0.15 (mol/mol), default_value=0.15, description="P:A molar ratio"
  • P_limit[Float64]=1.0e-6 (mol m-3), default_value=1.0e-6, description="limiting P concentration below which co-precipitation is inhibited"
  • P_components[Vector{String}]=["-1\P", "TAlk"], default_value=["-1\P", "TAlk"], description="totals or components that should be modified when P is consumed from solution"

Methods and Variables for default Parameters

  • do_PA_coprecip
    • rate_PA_coprecip (mol P yr-1), VT_ReactProperty, description="rate of P co-precipitation"
    • P_conc (mol m-3), VT_ReactDependency, description="P concentration"
    • A_formation_rate_1 (mol m-3 yr-1), VT_ReactDependency, description="substance A formation rate"
  • RateStoich_rate_PA_coprecip
    • rate_PA_coprecip (mol P yr-1), VT_ReactDependency, description="rate of P co-precipitation"
    • [P_sms] (mol yr-1), VT_ReactContributor, description="generated by RateStoich rate=rate_PA_coprecip"
    • [TAlk_sms] (mol yr-1), VT_ReactContributor, description="generated by RateStoich rate=rate_PA_coprecip"
    • [PA_sms] (mol yr-1), VT_ReactContributor, description="generated by RateStoich rate=rate_PA_coprecip"
  • totals
    • rate_PA_coprecip (mol P yr-1), VT_ReactDependency, description="rate of P co-precipitation"
    • rate_PA_coprecip_total (mol P yr-1), VT_ReactProperty, description="total rate of P co-precipitation"
source
PALEOaqchem.CoPrecip.ReactionPAReleaseType
ReactionPARelease

Release of P (eg iron-bound phosphorus) with A (eg Fe oxide) destruction

P=A -> P

at a rate Prelease = theta * A_destruction_rate_<n>, where theta = PA_conc / A_conc

P_components defines totals or components that should be modified when P is released (eg ["P"] to return P to solution).

Parameters

  • A_rate_stoich_factors[Vector{Float64}]=[1.0], default_value=[1.0], description="stoichiometry factor to multiply each A destruction rate variable to convert to mol A"
  • P_components[Vector{String}]=["P", "-1\TAlk"], default_value=["P", "-1\TAlk"], description="totals or components that should be modified when P is released"

Methods and Variables for default Parameters

  • do_PA_release
    • rate_PA_release (mol P yr-1), VT_ReactProperty, description="rate of coprecipitated P dissolution"
    • PA_conc (mol m-3), VT_ReactDependency, description="adsorbed P concentration"
    • A_conc (mol m-3), VT_ReactDependency, description="adsorbant concentration"
    • PA_theta (mol/mol), VT_ReactProperty, description="P:A molar ratio of adsorbed of coprecipitated P"
    • A_destruction_rate_1 (mol m-3 yr-1), VT_ReactDependency, description="substance A destruction rate"
  • RateStoich_rate_PA_release
    • rate_PA_release (mol P yr-1), VT_ReactDependency, description="rate of coprecipitated P dissolution"
    • [PA_sms] (mol yr-1), VT_ReactContributor, description="generated by RateStoich rate=rate_PA_release"
    • [P_sms] (mol yr-1), VT_ReactContributor, description="generated by RateStoich rate=rate_PA_release"
    • [TAlk_sms] (mol yr-1), VT_ReactContributor, description="generated by RateStoich rate=rate_PA_release"
  • totals
    • rate_PA_release (mol P yr-1), VT_ReactDependency, description="rate of coprecipitated P dissolution"
    • rate_PA_release_total (mol P yr-1), VT_ReactProperty, description="total rate of coprecipitated P dissolution"
source